library(terra)
library(tidyverse)
library(ggspatial)
library(sf)
library(patchwork)

# dictionary
#cc = current suitability
#ff = future suitability

# load shapefile
site = read_sf(file.choose()) %>% filter(District_N=='Mbinga')

site = read_sf(file.choose()) %>% filter(District_N=='Kyerwa')

#Load celevation raster
dem= rast(file.choose())

# Load current suitability raster
dat= rast(file.choose())

# convert raster into points
cc = terra::as.points(dat) %>% sf::st_as_sf()

cc1 = terra::extract(dat,cc)

cc1 = cc1[-c(1)]

names(cc1) = 'weight'

# convert sf points object to the dataframe
cc2 = cc %>% st_coordinates() %>% as.data.frame()

# merge pp2 and pp1
cc3 = cbind(cc2,cc1)

# calcuate for centroid
cweight_x_centroid = sum(cc3$X*cc3$weight)/sum(cc3$weight)
cweight_y_centroid = sum(cc3$Y*cc3$weight)/sum(cc3$weight)

cweighted_centroid = data.frame(x=cweight_x_centroid,y=cweight_y_centroid,period='Current')

ccentroid = st_as_sf(cweighted_centroid,coords=c('x','y'),crs=4326)

#----------------------------------------------------------
# 2050s
dat1 = rast(file.choose())

# convert raster into points
ff = terra::as.points(dat1) %>% sf::st_as_sf()

ff1 = terra::extract(dat1,ff)

ff1 = ff1[-c(1)]

names(ff1) = 'weight'

# convert sf points object to the dataframe
ff2 = ff %>% st_coordinates() %>% as.data.frame()

# merge two dataframes
ff3 = cbind(ff2,ff1)

# calcuate for centroid
fweight_x_centroid = sum(ff3$X*ff3$weight)/sum(ff3$weight)
fweight_y_centroid = sum(ff3$Y*ff3$weight)/sum(ff3$weight)

fweighted_centroid = data.frame(x=fweight_x_centroid,y=fweight_y_centroid,period='2050s')

fcentroid = st_as_sf(fweighted_centroid,coords=c('x','y'),crs=4326)

#-------------------------------------------------------------------------------

# 2070s
dat5 = rast(file.choose())

# convert raster into points
ff5 = terra::as.points(dat5) %>% sf::st_as_sf()

ff15 = terra::extract(dat5,ff5)

ff15 = ff15[-c(1)]

names(ff15) = 'weight'

# convert sf points object to the dataframe
ff25 = ff5 %>% st_coordinates() %>% as.data.frame()

# merge two dataframes
ff35 = cbind(ff25,ff15)

# calcuate for centroid
f5weight_x_centroid = sum(ff35$X*ff35$weight)/sum(ff35$weight)
f5weight_y_centroid = sum(ff35$Y*ff35$weight)/sum(ff35$weight)

f5weighted_centroid = data.frame(x=f5weight_x_centroid,y=f5weight_y_centroid,period='2070s')

f5centroid = st_as_sf(f5weighted_centroid,coords=c('x','y'),crs=4326)


#-------------------------------------------------------------------------------

# change of suitability change
dat2 = dat1-dat

dat21 = as.data.frame(dat2,xy=T)

dat3 = dat5-dat

dat31 = as.data.frame(dat3,xy=T)

# extract elevation values per points
dem_extract = terra::extract(dem,cc)

dem_extract = dem_extract[c(-1)]

names(dem_extract) = 'elevation'

# quantify change in suitability index between current and future along elevation gradient
c1 = cbind(cc1,dem_extract) %>% mutate(period='Current',scenario='Current')
f1 = cbind(ff1,dem_extract) %>% mutate(period='2050s',scenario='SSP245')
f2 = cbind(ff15,dem_extract) %>% mutate(period='2070s',scenario='SSP245')

dd = rbind(c1,f1,f2)

names(dd) = c('suitability','elevation','period','scenario')

dd = cbind(dd,cc2)

setwd('C:/Users/Mr. Kasongi/Documents/SDM paper revisions and new results')

write.csv(dd,'change in suitability along elevation ssp245 Robusta.csv')

# distance analysis of weighted centroid
dist1 = st_distance(ccentroid,fcentroid) %>% as.numeric() %>% as.data.frame() %>% mutate(period='2050s')
names(dist1) = c('distance','period')
dist2 = st_distance(ccentroid,f5centroid)%>% as.numeric() %>% as.data.frame() %>% mutate(period='2070s')
names(dist2) = c('distance','period')

distance = rbind(dist1,dist2)

setwd('C:/Users/Mr. Kasongi/Documents/SDM paper revisions and new results')

write.csv(distance,'change in suitability distance ssp245 Robusta.csv')

#ggplot(dd,aes(elevation,suitability,color=period))+geom_smooth()

#visualize results
p0 = ggplot()+geom_tile(data=dat21,aes(x=x,y=y,fill=mean))+
  geom_sf(data=site,fill='transparent')+
  geom_sf(data=ccentroid,aes(color=period))+geom_sf_text(data=ccentroid,aes(label=period),vjust=0.9,color='white')+
  geom_sf(data=fcentroid,color='red')+geom_sf_text(data=fcentroid,aes(label=period),vjust=-0.5,color='white')+
  geom_segment(aes(x=cweighted_centroid$x,y=cweighted_centroid$y,xend=fweighted_centroid$x,yend=fweighted_centroid$y),
               arrow = arrow(length = unit(0.2, "cm"),type = 'open'),
               size=0.5)+
  labs(x='',y='',fill='Change in coffee suitability scores',title = '(A)')+
  scale_fill_viridis_c()+
  scale_color_manual(values = 'red',label='Weighted centroid')+
  scale_x_continuous(breaks = seq(34.7,35.3,0.2))+
  theme_bw()+
  theme(legend.position = 'none')+
  annotation_north_arrow(location='tl',height = unit(0.7,'cm'),width=unit(0.7,'cm'))+
  annotation_scale(location='br')


p1 = ggplot()+geom_tile(data=dat31,aes(x=x,y=y,fill=mean))+
  geom_sf(data=site,fill='transparent')+
  geom_sf(data=ccentroid,aes(color=period))+geom_sf_text(data=ccentroid,aes(label=period),vjust=0.9,color='white')+
  geom_sf(data=f5centroid,color='red')+geom_sf_text(data=f5centroid,aes(label=period),vjust=-0.5,color='white')+
  geom_segment(aes(x=cweighted_centroid$x,y=cweighted_centroid$y,xend=f5weighted_centroid$x,yend=f5weighted_centroid$y),
               arrow = arrow(length = unit(0.2, "cm"),type = 'open'),
               size=0.5)+
  labs(x='',y='',fill='Change in coffee suitability scores',color='',title = '(B)')+
  scale_fill_viridis_c()+
  scale_color_manual(values = 'red',label='Weighted centroid')+
  scale_x_continuous(breaks = seq(34.7,35.3,0.2))+
  theme_bw()+
  annotation_north_arrow(location='tl',height = unit(0.7,'cm'),width=unit(0.7,'cm'))+
  annotation_scale(location='br')


p2 = p0+p1

p2

setwd('C:/Users/Mr. Kasongi/Documents/SDM paper revisions and new results')

ggsave('crop migration Robusta ssp585.png',p2,height = 6,width = 9,units = 'in')

# Arabica north arrow settings: height = unit(0.9,'cm'),width=unit(0.9,'cm')

# clear the environment and plotting panel
rm(list = ls())
dev.off()
